/*
    Copyright 2007-2011 Patrik Jonsson, sunrise@familjenjonsson.org. 
    Based on code by Volker Springel, used with permission.

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file 

    Declaration of the multiphase class, which calculates the
    parameters of the Springel & Hernquist 2003 multiphase model.
*/

// $Id$

#ifndef __multiphase__
#define __multiphase__

#include "cooling.h"

/** Implements the multiphase model from Springel & Hernquist 2003.
    Depending on density, a certain fraction of the gas is assumed to
    be in cold clouds while the rest is in a diffuse hot phase.  This
    class can calculate the mass in volume fraction of the cold phase
    as well as the temperature of the hot phase.  All quantities are
    in CGS units.  */
class multiphase {
private:
  /// Boltzmann's constant
  static const double k = 1.3806e-16;
  /// proton mass
  static const double m_p = 1.6726e-24;
  /// gamma of the gas
  static const double gamma = (5.0/3);
  /// hydrogen mass fraction
  static const double x_H = 0.76;

  /// Cooling object for calculating cooling times.
  cooling c;

  /// mean molecular weight for neutral gas.
  const double meanweight_n;
  /// mean molecular weight for fully ionized gas.
  const double meanweight_i;

  // multiphase model parameters
  const double T_c;
  const double T_SN;
  const double A0;
  const double beta;
  const double t0_star;

  // derived quantities
  const double u_c;
  const double u_SN;
  const double rho_th;

  /* Compute the density threshold. */
  double compute_threshold() const;

public:
  multiphase() :
    meanweight_n(4 / (1 + 3 * x_H)),
    meanweight_i(4 / (8 - 5 * (1 - x_H))),
    T_c(1000.0),
    T_SN(3e8),
    A0(3000),
    beta(0.1),
    t0_star(2.1e9*3.15e7),
    u_c(k * T_c / 
	(meanweight_n * (gamma-1) * m_p)),
    u_SN(k * T_SN / 
	 (meanweight_i * (gamma-1) * m_p)),
    rho_th(compute_threshold())
  {}

  /** This constructor explicitly specifies the timescale for star
      formation and the threshold density, so we can set them to be
      consistent with our simulations. The time is specified in years
      and the density in solar masses per cubic kpc. */
  multiphase(double ts, double rt) :
    meanweight_n(4 / (1 + 3 * x_H)),
    meanweight_i(4 / (8 - 5 * (1 - x_H))),
    T_c(1000.0),
    T_SN(3e8),
    A0(3000),
    beta(0.1),
    t0_star(ts*3.15e7),
    u_c(k * T_c / 
	(meanweight_n * (gamma-1) * m_p)),
    u_SN(k * T_SN / 
	 (meanweight_i * (gamma-1) * m_p)),
    rho_th(rt*6.77e-32)
  {}
    
  /* Returns the cold mass and volume fraction at the given density,
     given in cgs units. */
  std::pair<double, double> cold_fraction(double density) const;
  /* Returns the temperature of the hot phase gas for the specified
     density, given in cgs units. */
  double hotphase_temperature(double density) const;
};

#endif
